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Baryon distribution amplitudes (DAs) are crucial for the theory of hard exclusive reactions. We 
present a calculation of the first few moments of the leading-twist nucleon DA within lattice QCD. 
In addition we deal with the normalization of the next-to-leading (twist-four) DAs. The matrix 
elements determining the latter quantities are also responsible for proton decay in Grand Unified 
Theories. Our lattice evaluation makes use of gauge field configurations generated with two flavors 
of clover fermions. The relevant operators are renormalized nonperturbatively with the flnal results 
given in the MS scheme. We flnd that the deviation of the leading-twist nucleon DA from its 
asymptotic form is less pronounced than sometimes claimed in the literature. 



I. INTRODUCTION 



The notion of baryon distribution amplitudes (DAs) refers to the valence component of the Bethe-Salpeter wave 
function at small transverse separations and is central for the theory of hard exclusive reactions involving baryons [1, 
2, 3, 4, 5, 6, 7, 8, 9]. As usual for a field theory, extraction of the asymptotic behavior (in our case for vanishing 
transverse separation) introduces divergences that can be studied by the renormalization-group (RG) method. The 
distribution amplitude (p thus becomes a function of the three quark momentum fractions Xi and the scale that serves 
as a UV cutoff in the allowed transverse momenta. Solving the corresponding RG equations in leading logarithmic 
accuracy [10, 11] one is led to the expansion 

n=o ;=o V^slMoj/ 

The summation goes over all multiplicatively renormalizable operators built of three quarks and n derivatives and /3o 
is the first coefficient of the beta function. The polynomials Pni{xi) and anomalous dimensions 7„; are obtained by 
diagonalizing the mixing matrix for the three-quark operators 

{Dl'q){Dl'q){DX'q) , k, + + ks ^ n , 

and the c„/(/Lto) are the corresponding (nonperturbative) matrix elements. 

The theory of nucleon DAs has reached a certain degree of maturity. In particular the scale dependence is well 
understood [12, 13] and it reveals important symmetries of the quantum theory that are not seen at the level of the 
QCD Lagrangian [14]. At the same time, they are much less studied as compared to the usual parton distributions. 
One reason is that the approach to the perturbative factorization regime in hard reactions appears to be slow. There 
is overwhelming evidence that, e.g., electromagnetic and transition form factors at currently available momentum 
transfers of the order of a few GeV^ [15, 16, 17, 18] receive large nonfactorizable contributions from large transverse 
distances, usually referred to as soft (Feynman) or end-point contributions, and possibly from higher-twist corrections. 
This is indicated, for example, by the fact that the helicity selection rules are strongly violated. Another reason is 
that nucleon DAs enter physical observables in a rather complicated way through convolution integrals, integrated 
with smooth functions of the momentum fractions. This makes an experimental determination of the DAs pointwise 
in Xi very difficult. A qualitative picture suggested by QCD sum rule calculations is that the valence quark with 
the spin parallel to that of the proton carries most of its momentum [8, 9, 19]. It is timely to make this picture 
quantitative; lattice QCD ist best suited for this purpose [20, 21], allowing us to evaluate nonperturbative hadronic 
matrix elements of local operators that enter the expansion in (1) in a fully controllable fashion, at least in principle. 

In this work we report on the calculation of the first few moments of the leading-twist nucleon DA and also the 
normalization of the next-to-leading (twist-four) DAs [22] using two dynamical fiavors of clover fermions. The reason 
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why we also consider higher-twist DAs is that they enter the calculation of the helicity-violating Pauli form factor of 
the nucleon in perturbative QCD [^.;] and also the calculation of the soft (end-point) corrections to the form factors 
in the framework of the light-cone sum rule approach [24, 25]. Their knowledge is imperative for a QCD description 
of exclusive reactions in the JLAB energy range. It turns out that the same matrix elements are responsible for 
proton decay in Grand Unified Theories (GUTs), so they are also interesting in a broader physics context. A short 
presentation of our main results has already been given in Ref. [26, 27]. 

The paper is organized as follows. Section II contains a brief review of the general framework and definitions of the 
specific quantities that will be calculated. We focus on the relations to local matrix elements including those that are 
relevant for proton decay. 

In Section III we explain the lattice approach to the calculation of the matrix elements. The advantages of this 
method come at the cost of reduced symmetry due to the discretization of space-time. This leads to additional 
(unwanted) operator mixing as compared to the continuum, which has to be reduced as much as possible by a suitable 
choice of the operator basis. In particular, mixing with lower-dimensional operators is dangerous. The theoretical 
basis for the corresponding analysis is the classification of operators according to irreducible representations of the 
relevant lattice symmetry group. For quark-antiquark operators such a classification has been worked out in Ref. [28], 
while the analogous classification for the three-quark operators needed here is treated in Refs. [29, 30]. 

Section IV is devoted to the presentation of the numerical results for the matrix elements. We apply two different 
methods to analyze the data. The first one, which we refer to as unconstrained, is used to determine the normalization 
constants and to check the consistency of our results for higher moments. In the second method we use the momentum 
conservation as an additional constraint. This allows us to improve the accuracy of our results for the higher moments. 

In Section V we construct a model for the leading-twist DA, presenting our results in form of the canonical expansion 
Eq. (1), and compare it with other models in the literature. The final Section VI is reserved for a summary and 
conclusions. 

Some further technical details are presented in the Appendices, in particular the relations between the local operators 
relevant for leading-twist DAs of spin-1/2 baryons and the irreducible three-quark operators. We also present here 
the bare lattice results. 



II. GENERAL FRAMEWORK 



A. Leading twist 



The leading-twist proton DA can be defined [31, 32] from a matrix element of a gauge-invariant nonlocal three-quark 
operator: 



|u° {Zl)up {Z2)(f^ (2:3) Ua'a{zi, ZQ)Ub'b{z2, ZQ)Ub'h{zz, Zfi)e°' ''\p) 

= In 
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(2) 



Here cr^i/ = | [7^^, 7i/], C is the charge conjugation matrix, \p) is a proton state with momentum p, and N is the proton 
spinor; ellipses stand for the higher-twist constributions. All interquark separations are assumed to be light-like, e.g., 
u{zi) denotes the u-quark field at the space point zin with = 0, and U{zmZQ) denotes the non-Abelian phase 
factor (light-like Wilson line) 



U{zn,z,^) = P exp 



ig / dt [zn - zq) n^j^A^{tZn + (1 - t)zo) 



(3) 



Because of the light-cone kinematics, the matrix element does not depend on zq and the phase factors can be eliminated 
by choosing a suitable gauge. 

The invariant functions V , A and T can be presented in the form 

(4) 



(5) 



V{zip) = J[dx]exp - «^a;jZj(p • n) V{xi), 
and similarly for A and T, where the integration measure is defined as 

J [dx] = J dxi dx2 dx3 S{1 — xi — X2 — X3) . 
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The variables Xi have the meaning of the longitudinal momentum fractions carried by the three quarks in the proton, 
< Xi <1 and ^Xi = 1. 

The identity of the two u-quarks in (2) implies the following symmetry properties [32] 

y(xi,a;2,a;3) = y(a;2,a;i,X3), A(a;i, a;2, 2:3) = -A(x2, xi, X3), T{xi,X2,X3) ^ T{x2,xi,X3). (6) 

In addition, the requirement that the proton has isospin 1/2 yields the relation 

2T(xi,X2,X3) ^[V - A](xi,X3,X2) - [V - A]{x2,X3,Xi) (7) 

so that all three invariant functions can be expressed in terms of a single DA ip defined as 

(p{xi,X2,X3) = V{xi,X2,X3) - A{xi,X2, X3) . (8) 

The normalization convention is such that 

[dx]ip{xi,X2,X3) ^ 1 . (9) 



The definition in (2) is equivalent to the following form of the proton state [9, 32] 

b,T)=/~ / :^^^itL{\uHx,)uHx2)dHx3))-\uHx,)dH^2)uHx3))}, (10) 
J 2^J2^XxX2X3 

where the arrows indicate the helicities and the standard relativistic normalization for the states and Dirac spinors is 
implied. 

Moments of DAs are defined as 

y''""- C\Ax\x\x^xlV{x^,X2,X3) (11) 
Jo 

and similarly for the other functions. They can be related to matrix elements of the local operators 

yprmnj^Q-j ^yp(Ai---Ai)(/ii---/i,„)(iyi---i/„)|-Q-| 
^-p'nin^Qj ^,^p(Ai---A,)(aii---M™)(i'i---i'„)^q-j 

^e'^^^ {^D^- . . . D^^umi (C(-za«''))^^ {i^D^- . . . i?^'"7.(0)]^ {i^^D^- . . . i?''" (7^75^(0))]^ (14) 

by 

Pltw (0|V,^''^^(0)b) = -/ArF"">Vp>«iV,(p), (15) 

PlTW (0|^^''^«(0)b) = -/ArA""VpW"'A^r(p), (16) 

Pltw (0|r/''^«(0) \v) = 2/^^T"">Vp">«7V,(p) . (17) 

In the following we refer to these local operators as DA operators in order to distinguish them from three-quark 
operators with a general spinor index structure. The multi-index Ifhn with ^ = Ai . . . A/ (and similarly for m and n) 
denotes the Lorentz structure given by the covariant derivatives — dfi, — i%A^ on the right-hand side of Eqs. (12)- 
(14). The indices l,m,n (without bars) are the total number of derivatives acting on the first, second and third 
quark, respectively. A certain moment, e.g., y'™"^ is related to several operators which differ only by their 

Lorentz indices. Therefore the moments y'™"^ j^imn ^ rpimn ^.j^g right-hand side of Eqs. (15)-(17) can be calculated 
from different operators with same number of derivatives acting on the quark fields. The index p corresponds to the 
uncontracted Lorentz index of the gamma matrices in the operators. The leading-twist projection, Pltw^ can be 
achieved, e.g., by symmetrization in Lorentz indices and subtraction of traces. Our approach for handling the reduced 
symmetry of the discretized space-time properly is described in Section III. 
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The symmetry relations (6) are translated into similar relations for the moments: 

-ylmn -^m/n ^/mn j^mln rj-ilmri rjiirdn (18) 

For further use we define the combination 

^Imn ^ ^-yhnn ^Imn _|_ Q^rj-ilnrrfj (19) 

3 

Taking into account the isospin relation (7), the moments of V, A, T can be restored from the moments of cj) by 

T'"" =^(0'"™ + 0™"'), (20) 

ylmn ^1 (20'™" + 2(/>"'" - 0"'" - 0""') , (21) 

A'"" =^ (-20'"" + 20"'" - 0"'™ + 0""') . (22) 
The conventional proton DA (p{xi) (8) is given in terms of 0(a;i) as 

V{xi,X2,X3)^2<j>{xi,X2,X3)-cj){x3,X2,Xi), (p'™" = 20'"" - 0""' . (23) 

Due to momentum conservation {xi + X2 + X3 = 1) there are additional relations between lower and higher moments: 

0'"™ = _|_ _|_ 0''n("+l)^ (24) 

In particular this implies 

1 = 000" = 0100 + 0010 + 0001 ^ ^200 ^ 0020 ^ ^002 ^ 2(00" + ^101 + ^hO) ^ (25) 

B. Next-to-leading twist operators and proton decay 

In general, there exist three independent next-to-leading (twist-four) three-quark DAs, cf. Ref. [22]. In this work we 
only consider their normalization, which is related to the contributions of local operators without derivatives. Thus 
the problem is simplified considerably since the general Lorentz decomposition of the relevant matrix element involves 
only four structures: 

4(0|e"''X(0)^^(0)rf;;(0)b) ^V,%fCUi^,N)^ + V^mNh^CUiri^N^ 

+ T0(p^ia^,C)„^(7^75^)^ + T0m^(a^,C)a/3(a''"75iV)^ , 

where uin is the nucleon mass and we have used the same notation as in [22]. The leading-twist-three constants Vi 
and Ti correspond to /atV^"'^*' and /atTOOO -^^ q^j. notation, Eqs. (15) and (17), and are equal. The two additional 
constants, Vg" and Tg , correspond to subleading twist-four contributions. The combinations Ai = — iV^ and 
A2 = Q{Vi — 4T0) are often arising in QCD sum rule calculations. They describe the nucleon coupling to the two 
independent local operators 

£,(0) = e-"' [u''^{0)CYu\0)] X {j^jpd%0))r, (27) 

Mr{0) = e-'^ [u''^{0)Ca''''u\0)] x (75(7^.^^(0)), , (28) 

which have been introduced in [33, 34], respectively. Their matrix elements are given by 

{0\Crmp) ^ Xim^Nr, (29) 
{0\MrmP) ^ >^2mNNr. (30) 

Separating the components of different helicity, one can write 

Cr=4{jBU'^ -JlU^^X , (31) 

Mr^ShnU'^-jLU'^). (32) 
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where 7l = (1 75)/2 , 7^ = (1 + 75)/2 are the left- and right-handed projectors and 



^abca 



ilL/Ru'') C-fL/Rd' 



The Fierz identity impHes 

-"^^ [u''^{0)Cru\0)] (757m'^'(0)). = 2e'"^ (- [u''^ {0)Cj,d\0)\ u-{0)r + [u°^(0)Cd''(0)] (75"=(0)), 



e 

^abc 



'u''^{0)Ca^^''u''{0)\ (75^^.rf^(0)), = 4e''^^ ([w'^^(0)C75d'(0)] u^(0), [u''^ iO)Cd\0)\ (75^^^(0)), 
Thus we get 

mAr(2Ai + A2)iV(p) = 8(0|e'^''^ (u^^Cd") 75"1p) , 



(33) 

(34) 
(35) 

(36) 



where, as it can be shown, the matrix element on the right-hand side vanishes in the nonrelativistic limit. 

The operators (27) and (28) appear also in the low-energy effective action of generic GUT models, and their matrix 
elements {Tr\£\p) and (7r|A^|p) give rise to proton decay. These matrix elements, in turn, can be related to the 
constants defined in (29), (30), using soft pion theorems or, what is the same, leading order in chiral perturbation 
theory [35, 36, 37, 38, 39]. 
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FIG. 1: Diagrams contributing to the nucleon decay amplitude p ^ tt" + e"'". 

To this end one introduces two low-energy constants a and /? which extend the usual three-flavor baryon chiral 
Lagrangian. They are defined by 



(0|(7LZ^^).(0)b) = -a(7LiV),, 
(0|(7LZ^^).(0)b) = -/3(7LA^)r, 



(0|(7flZ^^).(0)b) = a(7i?iV),, 
(0|(7flW^).(0)b)=/3(7fliV),. 



(37) 
(38) 



Because of (31), (32) one obtains a = itinXi/A and f3 = 171^X2/8. The knowledge of these two constants allows one 
to estimate nucleon-to-pion decay matrix elements. Using the notation of Ref. [40] the relevant factors in the decay 
amplitude for the proton to 7r° decay (cf. Fig. 1) have the form 



<^(p-0 



V2f 



(1 + ffA), 
(1 + ffA), 



(39) 
(40) 



where / is the tree level pion decay constant normalized such that the experimental value is /^r ~ 131 MeV and gA is 
the axial charge. 



III. DETAILS OF THE LATTICE CALCULATION 



In this section we discuss the techniques used and the details of the lattice calculation. From now on we work in 
Euclidean space. In order to define the Euclidean counterparts of the operators Eqs. (12)-(14), (27)-(28) and (33) we 
interpret the Dirac matrices and coordinates as being Euclidean. For our Euclidean Dirac matrices see Appendix A. 
The expressions on the right-hand-side of Eqs. (12)-(14), (27)-(28) are then modified accordingly. In the first part of 
this Section we summarize the general features of our approach. The following parts contain the description of the 
calculation of matrix elements relevant for leading and next-to-leading twist DAs. 
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A. General features 

To be as flexible as possible in our calculation we have adopted a two-stage approach in the evaluation of the 
correlators. In the first step we have calculated correlators of the form 

CTp^^ = {e^'^[D^,...D>^M^)]l[D,,...D,^u{x)%[D,,...D^^^^^^ (41) 
with / + m + n < 2. As interpolating operator for the proton we have used 



Mr 



u'^^C-fsd'' u%. (42) 



Due to the presence of two u-quarks in the three-quark operator, C^"^ can be reconstructed from C^™"t- by an 
appropriate interchange of Dirac indices. 

In the second step the general three-quark operator from Eq. (41) was used to calculate the matrix elements for 
the different quantities we discussed before. The general form of the correlation functions we compute at this stage 
reads after projection onto momentum p: 

{Or{t,p}Afr'{0,p)) - J2^0\Ormp,s)Nr,{p,s)exp{-E{p)t)- (43) 

Here contributions of excited states have been neglected and the dependence of the nucleon states and spinors on the 
spin vector s has been made explicit. For the energy E{p) we use the continuum expression E{p) = \/m% + p'^. We 
have checked that this dispersion relation is fulfilled well within errors (see, e.g.. Fig. 2), so we had to fit only the 
mass in the exponential. The correlator in Eq. (43) can be directly constructed from the general correlation function 
(41). The matrix element on the right-hand side is the quantity we want to determine. Thus we have also to calculate 
the normalization constant Zn{p), which can be extracted from the usual two-point nucleon correlator 

Cn{p) ^ (7+).v {Afrit, mr'{0,P)) = ^^(P) "''^j;f^^^^ exp {~E{p)t) (44) 

E{p} 

with the positive parity projection 7+ = (1 -I- 74)/2. In the evaluation of the correlator in Eq. (41) the overlap of 
the nucleon interpolator with the nucleon state is improved by Jacobi smearing at the source while the sink is not 
smeared since we want to evaluate local matrix elements. Thus the nucleon correlator in Eq. (44) cannot be extracted 
from the general three-quark nucleon correlator (41) but must be computed separately with Jacobi smeared sink and 
source. 

The normalization constant Zjv(p) could be removed by considering the ratio 

{i^+)r'r{Ormr^mf 

(7+),,,(AA,(t)Ar,,(o)) ■ ^ ' 

However, as we will see later, the location of the effective mass plateaus is different for the two correlators, presumably 
due to the different smearings on the sink, spoiling this simple approach. Thus instead of calculating the ratio we 
perform a correlated fit to the two correlators in the range of the corresponding effective mass plateaus. 

Up to now we did not take into account that our calculations are performed on a space-time lattice. This leads 
to reduced symmetry compared to the continuum. Due to this symmetry reduction we expect additional operator 
mixings which are not present in the continuum. In particular, we can have mixing with lower-dimensional operators. 
Thus a systematic analysis and careful choice of the operators used is mandatory. In [30] a complete classification 
with respect to the spinorial extension of the hypercubic group H{A) for all three-quark operators without derivatives 
is presented. For operators with one and two derivatives the classification is worked out for the leading-twist case. 
These results enable us to derive operators with "good" mixing properties, good in the sense that they do not mix 
with lower-dimensional operators. They belong to definite irreducible representations of the spinorial extension of 
i/(4) and are most easily constructed in the Weyl representation of the Dirac matrices. Therefore we also work in 
this representation. 

In Table I we give an overview of the irreducible multiplets of operators taken from Table 4.1 in [30], with a modified 

notation adapted to our needs, e.g., operator B^J corresponds to cj^^^j^ in [30] and similarly for the others. The 

next-to-leading twist operators (27) and (28) lie completely within the Tf representation with mass dimension 9/2. 
The operators relevant for the leading-twist DAs belong to other multiplets. As operators without derivatives in the 
T- representation do not have an overlap with the nucleon, the relevant operators with "good" mixing properties lie 
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d = 9/2 (0 derivatives) 


d = 11/2 (1 derivative) 


d = 13/2 (2 derivatives) 


4 


r{0) k{0) „(0) „{0) „{0) 
^1 i 1 ^2 i 1 ^3 i ' ^4 i ' ^5 i 




K(2) k(2) „(2) 
"l i ' "2 i 1 "3 i 


4 






k(2) r(2) „(2) 
"4,1 ' "5,1 ' "6,1 




g{0) 
"6,1 


g(i) 


r(2) »(2) „(2) 
"7,1 ' "8,1 ' "9,1 


12 


r{0) c{0) »{0) 
"7,i ' "8,i ' "9,i 


r(i) 

"2,i ' "3,1 ' "4,8 


r(2) c(2) c(2) »(2) 
"lO,i' "ll,i> "l2,i> "l3,i 


12 
"^2" 




"5,i ' "6,i ' "7,i 1 "8,1 


k(2) k;(2) „(2) „(2) „(2) 
"l4,ii "l5,i' "l6,i' "l7,i' "l8,i 



TABLE I: Overview of irreducibly transforming multiplets of three-quark operators sorted by their mass dimension (number 
of derivatives) taken from [30] with a notation adapted to our needs. Since for the classification it is not important on which 
quarks the derivatives act, only the sum I + rn + n is given as a superscript. The subscript gives the numbering of the operators 
according to the numbering convention in [ j ,]. The first number corresponds to the lower index of [.']()] while the second number 
corresponds to the upper index in [30] labelling different operators within one multiplet (cf. Table 4.1 in [30]). In the first 
column we give also the representations in the notation of [30] where the superscript denotes the dimension. 



in T^, and for zero, one and two derivatives, respectively. Rewriting these irreducible operators in terms of 
the DA operators defined in (12)-(14) allows us to choose those that are suited for lattice calculations. The ensuing 
relations for leading-twist spin-1/2 baryon DAs are summarized in Appendix B. In the following we give some details 
for these operators. 

Initially, the irreducible operators in [ ■ ] have a general flavor content. Considering the case of two derivatives as 
an example we have operators of the type 

Tf^-^D^D^e'^'^r^g^h'^^ , (46) 

where T'^^^'^ is a tensor projecting the operator to a certain irreducible representation. As it is not important for the 
construction of irreducibly transforming operators on which of the quarks the derivatives act, the different possibilities 
fall into the same irreducible representation. The proton operators are then recovered by the identification 

/ ^ li, g-^u, d, (47) 

and subsequent projection onto isospin 1/2, which is done by combining properly different multiplets. This procedure 
differs somewhat from the approach adopted in Ref. [30], but it leads to equivalent results. 

The operators used in our calculation have to be renormalized. In [41, 42] the required renormalization matrices were 
calculated nonperturbatively on the lattice imposing an RI'-MOM-like renormalization condition. Using continuum 
perturbation theory and the renormalization group the results were converted to the MS scheme at a scale of 4GeV^. 
Note that in this procedure the mixing with "total derivatives" is automatically taken into account. The scale at 
which our renormalization condition is imposed is taken to be 20GeV^, and the systematic uncertainty is estimated 
by varying this scale between lOGeV^ and 40GeV^. 



B. Moments of the leading-twist DA 



0th moment 



12 

Using the representation and the relations to the DA operators given in Appendix B we construct three 
quadruplets of operators with isospin 1/2 from the twelve irreducible three-quark operators, which can be used to 
calculate /at: 



( 



KjOOO 
""8,6 
rOOO 
A-'a 1 



-B\ 



odo 



12 
rOOO 



KjOOO 

"9.6 

rOOO 

"9.1 

rOOO 

"9,12 

rOOO 

"9,7 







+ Br, \ 


4 




rOOO 
"9.3 


3 


Kjodo 

"8.10 


+ BtZ 




I K7000 

\ "8,9 


- S9°°9° ) 



4V5 



"8,5 
rOOO 
"8,8 

\ nom 

\ "8,11 



BO™ ^^ 

"9,5 
"9,8 
"9,11 / 



(48) 
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The three-quark operators O on the left-hand side have also a Dirac index which we do not give explicitly here. The 
relations to the DA operators given in Appendix B yield then 



{Q\0T^\P) = fNiipili - ^P2l2)N{p), 

(0|O^™|p) = fNiipili + «P272 - «P373 + E{p)-(i)N{p). 



(49) 
(50) 
(51) 



The operators O^q and are most suitable for our calculation since '^ould require nonzero spatial momenta 
in the 1 or 2 direction, which would increase the statistical noise. Thus, in order to determine /jv, we evaluate finally 
only the following two correlators at p = (J: 



EE {{^,OT,{t,P))^ {m,P))r, {l+)r'r) 



pI 



expi-Eip)t): 



(52) 



C 



■000 
C,0 



fN^Z^ip) — 

fNVZN[p) exp{-E{p}t). (53) 



1st moments 



We use the irreducible operators with one derivative from Appendix B to construct operators for the calculation of 
the first moments of the proton DA, 



Iran 



//-j/mn 



A^2 



( Ki 



Qlnra \ 



6,2 
Imn 



f2t7m 

. f2lmn 
\ "6,8 



"7,7 
'->7 s 



A^/2 



/ K?lmn 
I ^6,3 



f2lmn 

Irnn 



f2lnm \ 
^7,3 \ 
f2lnm 

f2lnm 
^7,9 
f2lnm , 
^7,10 / 



/ 


K?/mn 
^6,6 


Y2lnm \ 
^7,6 \ 




rAmn 
^6,5 


f2lnm 
~ ^7,5 




r2lmn 
^6,12 


_i_ f2lnm 
~^ ^7,12 


V 


rAmn 
~^6,11 


+ ^7?n / 



(54) 



where the the superscript Imn with I 



1 and nonnegative integers l,m,n indicates on which fields the 



derivative acts. The matrix elements of these operators are then 

(0|O^™"b) [(pi7i -P272)(*P373 - E{p)j,) - 2zpip27i72] Nip), 

(0|O^":rb> =/iv0""" [{pUi +P272)(*P373 + E{p)j,) - 2p^E{p)j3l4] N{p), 
(0|O^™i"b) =/^0"""(-pi7i +P212){^P313 + E{p}j,)N{p), 



(55) 
(56) 
(57) 



where again a Dirac index is implied for the three-quark operators O. Unlike the case of the 0th moment all operators 
require at least one nonzero component of the spatial momentum. Hence using all operators available in this case we 
evaluate the correlators 



. /^T^ Eip^{mN + E{p^)+2pl~pl 

VZn[P^ Pi ^^7^^ cxp {-E{p)t) , 

E{p) 

(58) 



^ ((747iObT(^,P)). {m^,P))^, h+)r'r) = 
^ ((747lO^T(i,P)), {m,P))r, h+)r'r) - 

to determine the first moments (^^™, S'^-^'^ and 6'^'^^. 



fN^'^^^^/Z^p^ 
fNCP^'^'^VZ^Pl 



E{p}{mN + E{p}) + pI 

E{p) 
E{p){mN + E{p)) 



pI 



E{p) 



exp{-E{p)t), (59) 



exp {-E{p)t) (60) 



2nd moments 



The calculation of the second moments requires the use of the four-dimensional irreducible representation to avoid 
mixing with lower-dimensional operators. Unfortunately, this decreases also the number of possible operators. Using 
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the irreducible three-quark operators with two derivatives and the relations to the DA operators from Appendix B 
we construct 



3V3 



Inm 
6,4 


rAmn 
^5,4 


\ 


Inm 


f2lmn 




6,3 




Inm 
6,2 


- ^5:r 




Inm 
6,1 


f2lmn 


) 



where now / 



n + n — 2 with l,m,n nonnegative integers. The corresponding matrix element is given by 

(0|O^"|p) = [piP27i72 + E{p}ji) + ip3E{p}j3ji (zpi7i - N{p) 



and the second moments are determined from 



C^" EE ((7273740r"(i,p)), (Ar(0,p))^, (7+).,J = -/a.0"""v/^^ 



P2P3- 



E{p){mN + E{p)) 



p\ 



(61) 



(62) 



E{p) 



exp{~E{p)t)- 
(63) 



C. Next-to-leading twist DAs 



For the higher-twist DAs we consider only the operators without derivatives. If we write the operators in Eqs. (27)- 
(28) and (33) with general flavor content, 







Mr{0) 









r^{Q)C-ipg\Q)\ X (757p/i'(0)). , 
r^{Q)Ca^,g\Q)\ X (75^,../i=(0))r 

we can express them in terms of the irreducible three-quark operators as 

\ 



^Imn 

j^lmn 

Imn 
3,4 



"4,2 

f2lmn 

"4,3 

■ ^"1" 



M 



( Birr \ 

BiJT 

, f2lmn 
\ "2,4 



and 



7hZ^^ - jlU' 



After the identification f ^ u, g 



( 

^3^" 
f2lmn 
"3,3 

\ "3,4 



7iZ^^ - 7flZ^^ 



^273 



hn'fi 
1,1 


~ ^2,1 


\ 


hn7i 
1,2 


r2lm7i 
~ ^2,2 




^mn 
1,3 


rshnn 
~ ^2,3 




lm7i 

1,4 


O2 4 


/ 



h ^ d, we restore the proton operators in (27)- (28) and (33). 



(64) 
(65) 
(66) 



(67) 



(68) 



IV. NUMERICAL RESULTS 



We have evaluated our correlators on the QCDSF/DIK configurations generated with two flavors of clover fermions 
at two different (3 values summarized in Table II. For /3 — 5.29 we have used two different lattice sizes, 24^ x 48 and 
16"^ X 32, each at three different quark masses. For /3 = 5.40 we have evaluated the correlators at five different quark 
masses on 24^^ x 48 lattices. The lattice spacing has been set via the Sommer parameter rg = 0.467fm [43, 44]. As far 
as possible we have also checked that the dependence of the final results on the fitting procedures discussed below is 
only very mild and the deviations are consistent with the present statistical errors. 



A. General discussion 



As already anticipated we can reduce the noise by combining different momenta and/or different operators. However, 
calculating the general three-quark operator for many momenta turned out to be too expensive. Hence the general 
correlators (41), and therefore also the correlators for DA operators, were evaluated only for a minimal set of momenta. 
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p 




// ^TT [VjrC V J 


T Tf^ In TYl d 

VUl Llllic 


a[fm] 


L[fm] 


5.29 


0.1340, 0.1350, 0.1359 


1.411, 1.029, 0.587 


16^ X 32 


0.08 


1.28 


5.29 


0.1355, 0.1359, 0.1362 


0.800, 0.587, 0.383 


24^ X 48 


0.08 


1.92 


5.40 


0.135, 1356, 0.1361, 
0.13625, 13640 


1.183, 0.856, 0.648, 
0.559, 0.421 


24^ X 48 


0.07 


1.68 



TABLE II: The set of lattices used in our calculation. The scale was set via the Sommer parameter ro — 0.467 fm. 

To extract the nucleon wave function normalization constant we have fitted the correlator 

Cr-liC^o + Cro) (69) 
where we have averaged over the two possible correlators at p = 0. Similarly, for the first moments we have used 

1 — 2 V-'A,! + + '^C,l j I \'^) 

with I + m + n — 1 and p— (27r/L, 0, 0), where L is the spatial extent of our lattice. For the second moment we have 
only one correlator, hence no averaging is possible and we have evaluated it for p — (0, 27r/i, 2n/L). 

To determine the normalization constant Zm{p) we had also to evaluate the usual nucleon correlator. As the 
additional smearing on the sink introduces additional noise, in particular for p 0, we have improved the signal by 
using different momenta in the nucleon correlator. For the 16^ x 32 lattices we have worked with 

'^N = l {CN{2nlL, 0, 0) + Cjv(0, 2ti/L, 0) + Cn{Q, 0, 2n/L)) (71) 

and 

^N = \ {Cn{Q, 27r/L, + Cn{2tt/L, 0, 2n/L) + CAr(27r/L, 2tt/L, 0)) , (72) 

while for the 24^^ x 48 lattices we have used a larger number of momenta: 

'^N^\ {Cn{2^/L, 0, 0) + Cn{Q, 2n/L, 0) + Cjv(0, 0, 2^/L)) , (73) 

d =^ (CAr(0, 27r/i, 2tt/L) + (7^(0, -27r/i, 2tt/L) + Cn{2tt/L, 0, 27r/L) 

+ Cn{2t:/L, 0, -27r/L) + Cn{2t:/L, 2t:/L, 0) + Cn{2ttIL, -2t:/L, 0)) . (74) 

As already mentioned, the location of the effective mass plateaus for the nucleon correlator differs from that for the 
other correlators as exemplified in Fig. 2. Thus, instead of calculating the ratios of the correlators we have performed a 
joint fit. As all correlators are evaluated on the same gauge configuration we should also take into account all possible 
statistical correlations. We have employed two different fitting procedures with different possibilities for incorporating 
the correlations: 

PC: The first possibility is to fit every moment of the DA separately, e.g., for fN4>^^^ we fit the correlators Cl'^^ and 
C^f simultaneously and incorporate the correlations of both correlators and those between different time-slices. 
However, since we want to extract (/)^°° and not fN4>^'^^ we should in principle also consider the correlation with 
(^000 j-j^g ^jjg omission of these additional correlations we call this procedure "Partially Correlated". 

FC: For the second possibility we have estimated the full crosscorrelation matrix and call this method therefore 
"Fully Correlated". In this case we fit simultaneously the correlators for the zeroth, first and second moment as 
well as the nucleon correlator with the same modulus of the momentum. 

Both methods have some common disadvantages. In order to extract the moments we have to perform multiparam- 
eter fits which involve nucleon mass, different normalization constants and the moments. The second disadvantage 
is the required knowledge of the smeared-smeared nucleon correlator for nonzero spatial momenta, which introduces 
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.0,8 -J 
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(a) 
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PIG. 2: Effective energy plots for different nucleon momenta at /3 = 5.40 and k = 0.13610 for the nucleon correlator (a) 
and the distribution amplitude correlators (b), where we have averaged over all available correlators. The black circles were 
obtained at zero nucleon momentum, the red squares and blue diamonds correspond to p2 = (27r/L)2 and = 2{2tt/LY, 
respectively. The lowest black line shows the nucleon mass as obtained by direct calculation. The middle red and the top blue 
line correspond to energies iJ^H = Jti^h with = {2-k/LY and p^ = 2(27r/L)^, respectively. 
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FIG. 3: Plateaus of correlator ratios i?^"" (black diamonds) and ii^"" (blue squares) for /3 = 5.40 and k = 0.1361 together 
with the corresponding fit values and the associated error bands. 



additional noise. This requirement can be avoided if we consider ratios of the correlators, which are equal to ratios of 
moments: 



l+m + n = 2 



R 



Inin 



Si Css 



Si = 



/lOO I iOlO 



^ ^100 I ^010 I /-.OOl 



(75) 



So 



C. 



S,2 



S2 
Cs,2 



xOll I i,101 



/002 



2(C°" + 1 + C2"") + Cfo + Cf" + 2. (76) 



Now we need additional input to determine the normalization of the moments (p'-™-" with I + m + n > 1. This can be 
obtained by using the constraint (25). Thus, we require, e.g., for the first moments that the renormalized moments 
satisfy 



(77) 



12 




0,9 
ml [GcV 

(a) 



FIG. 4: Linear chiral extrapolation for Jn /m^ (a) and Xi/ruN (b) with the 1,2 and 3 sigma error bands. 

where (^J^' are the unrenormaUzed lattice values 

*,ioo 



and Z is the renormalization matrix. This leads immediately to a constraint for the ratios i?'^* — 

1 



3 ^1 



lat 



E'7 plat ■ 



(78) 



(79) 



As in this case we use explicitly the constraint (25) we call this analysis method "constrained". The calculation of the 
ratios i?'™" does not suffer from the disadvantages mentioned above. Fitting these ratios to a constant we can reach 
a much higher precision compared to the unconstrained method discussed before. In Fig. 3 we present some of these 
ratios obtained on one of the ensembles with /? = 5.40. They exhibit longer and less noisy plateaus compared to the 
correlators in Fig. 2. 

The lattice results are obtained at nonphysical quark masses and we have to extrapolate them to the physical 
point. To our knowledge there are no calculations in chiral perturbation theory to guide our extrapolation. Therefore 
we have to rely on the behavior of our data and extrapolate them linearly to the physical point. To estimate the 
systematic uncertainty of this chiral extrapolation we have performed also an extrapolation including a quadratic 
term. The systematic uncertainty is then taken to be the difference of the two results. 

In the following we present the results of the constrained and unconstrained analysis methods discussed before 
in the MS scheme at 4GeV^ while the raw lattice results are summarized in Appendix C. Using the unconstrained 
analysis we obtain the normalization constants of the DAs and test how good the constraint in (24) is satisfied. Better 
results with smaller errors for the higher moments of the leading-twist DA are then obtained from the constrained 
analysis. 



B. Unconstrained analysis 

In Table III we present the results for the different constants which are associated with operators without derivatives: 
the nucleon wave function normalization constant and the next-to-leading twist normalization constants Ai and A2. 
Our results confirm the relative signs of fjy, Ai and A2 calculated in ['}'}, 45]. Furthermore we observe TO7v(2Ai-|-A2) /8 = 
q; -I- /? « as in [40, 46]. This is expected since due to (36) 2Ai -I- A2 vanishes in the nonrelativistic limit and is known 
to be small at small quark masses [47]. 

Our results for the nucleon wave function normalization constant exhibit a clearly nonlinear behavior as a 
function of mj. However, the dimensionless ratio /Ar/m^ is approximately linear (see Fig. 4(a)) and it has the 
additional advantage that it does not suffer from the uncertainty in setting the scale on the lattice. The chiral 
behavior of Ai and A2 is less clear and we have performed two different chiral extrapolations for these quantities. First 
we have extrapolated the constants Ai linearly to the physical point and then we have applied the same procedure to 
the ratios Xi/niN. The linear fit looks more favorable for the ratios Xi/rriN (see Fig. 4(b)). Thus we take the results 
from this fit as our final values, but for comparison we also give the results from the other extrapolation. In contrast 
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to [40, 46] we do not observe linear behavior for rriNK as a function the quark mass. However, our results from the 
linear extrapolation of Xi/niN are compatible within the errors with those in [40, 46]. 

We have determined the moment combinations (^'™" — 2(/)'""' — 0"™' also directly and not from the results for 0'™", 
using the PC fitting procedure. Thus we had also to compute /jv within this approach. We have also determined 
using this analysis method. The results are presented in Table IV. The correlators for higher moments entering the 
FC fitting procedure seem to favor slightly larger nucleon masses, while the PC analysis leads to somewhat higher 
values of the normalization constants. We consider the values for the normalization constants obtained within the 
PC analysis to be more reliable as they are not perturbed by the noisier correlators for the higher moments. 



/3 


5.40 


5.29 


fN/m% ■ 10^ 


3.486(60)(56)(60) 


3.290(62)(100)(72) 


-Ai/jrijv • 10^[GcV] 40.64(65)(194)(110) 


41.24(72)(200)(128) 


-Ai • 10^[GoV^] 


49.84(95)(290)(135) 


52.47(104)(135)(164) 


Aa/mjv ■ 10^[GeV] 


80.17(131)(396)(218) 


82.08(146)(452)(254) 


A2 • 10^[GeV^] 


98.53(189)(601)(268) 


105.12(209)(250)(324) 




0.3457(75)(89)(3) 


0.3530(62)(132)(7) 




0.3124(81)(128)(4) 


0.3176(62)(108)(2) 




0.3142(77)(100)(4) 


0.3283(62)(68)(4) 


<^°" 


0.0838(73) (266) (44) 


0.0851(61)(1)(44) 




0.1121(92)(250)(58) 


0.1020(66)(179)(68) 




0.1051(67)(6)(4) 


0.0979(54)(5)(9) 


^200 


0.1523(106)(699)(129) 


0.1639(86)(216)(114) 


^020 


0.1268(97)(153)(98) 


0.1277(79)(1)(76) 


^002 


0.1398(99)(45)(128) 


0.1473(84)(40)(111) 



TABLE III: Chirally extrapolated results from the FC analysis for normalization constants and the moments 0'™" at /3 = 5.40 
and (3 = 5.29 in the MS renormalization scheme at 4GeV^. The first error is the combined statistical error of the moments 
and renormalization matrices. The second (third) errors are the systematic uncertainties due to the chiral extrapolation 
(renormalization) . 



f} 5.40 5.29 



/jv/mw • 


10^ 


3.672(78)(90)(63) 


3.538(79)(283)(77) 


-Ai/mjv ■ 10^[GeV] 42.19(81)(86)(115) 


45.07(92)(315)(140) 


X2/mN ■ 10^[GeV] 


82.91(171)(18)(225) 


86.90(87)(641)(261) 






0.3871(313)(528)(4) 


0.3903(204)(464)(12) 


= i 


^010 


0.3150(226)(290)(720) 


0.3298(159)(118)(608) 


001 




0.3155(272)(453)(2) 


0.3277(190)(270)(5) 






0.0712(180)(127)(92) 


0.0827(137)(103)(92) 




^101 


0.1091(112)(138)(64) 


0.1176(105)(171)(64) 






0.1266(178)(82)(40) 


0.1069(137)(103)(49) 






0.1879(250)(942)(135) 


0.1709(184)(569)(121) 


020 

V =4 


^020 


0.1275(149)(105)(108) 


0.1261(117)(78)(75) 


002 

V 




0.1357(233)(375)(135) 


0.1249(193)(296)(109) 



TABLE IV: Chirally extrapolated PC results for normalization constants and the moments i/j''"" at /3 = 5.40 and (3 — 5.29 in 
the MS renormalization scheme at 4GeV^. The first error is the combined statistical error of the moments and renormalization 
matrices. The second (third) errors are the systematic uncertainties due to the chiral extrapolation (renormalization). Note 
that only the values for (ys'™' can be directly compared with the values for 0'™' in Table III. 

As expected, the nonzero spatial momenta make the results for the first moments noisier than for operators without 
derivatives. The renormalized results for the moments 0^"°, (jp^^ and (^°°^ show clearly the deviation from the asymp- 
totic case with 0^"" = (fP^^ = (/i^*^^ = 1/3. As the relative differences of these moments describe the deviation from 
the symmetric case, they are of particular interest in phenomenological applications. Thus we have also determined 
these differences directly and the bare results from the PC analysis are given in Appendix C. Although these results 
show a significant deviation from the symmetric case, the errors are large and do not allow reasonable quantitative 
conclusions. To illustrate these we show in Fig. 5(a) the most important asymmetry 0^"° — (fp^^ normalized by the 
sum so that we can compare this later directly with the results from the constrained analysis. However, the results 
for the moments are less affected by the noise as shown on the example of (f)^^^ in Fig. 5(b) also normalized by Si. 

We have checked our results by calculating the sums Si and 5*2 according to Eqs. (75) and (76). The results for the 
bare and renormalized sums are shown Fig. 6. For the renormalized moments the constraint (24) is fulfilled very well 
indicating the consistency of our results. Of course the statistical and systematic errors for the case of two derivatives 
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FIG. 5: Chiral extrapolation of the asymmetry 



'-0°^°)/5'i (a) from PC results and the ratio ^^''V'S'i (b) from FC results. 
We have normalized the values by 5*1 so that we are able to compare these directly with the plots in the constrained analysis 
in Fig. 7. 



20 30 
^2 [GeV^] 

(a) 




FIG. 6: The bare (solid black line with statistical error band) and renormalized (blue diamonds) sum of the first moments 
(a) and second moments (b) according to Eq. (24) as obtained from the FC analysis. The smaller errors for the renormalized 
values are purely statistical, while the larger are ones include the systematical error due to the chiral extrapolation. The three 
different points were obtained from three different renormalization scales jj, in the RI'-MOM scheme to estimate the systematic 
uncertainty due to the renormalization. The theoretical constraint (25) that the sum should be exactly equal to one is fulfilled 
in both cases. 



in the operators are higher. Nevertheless, the results still allow us to see the asymmetries. Because of the large errors 
we give these only for the bare results in Appendix C. 



C. Constrained analysis of higher moments 



In the last section we have seen that the unconstrained analysis of our data gives us results consistent with theoretical 
constraint (25). However, better estimates of moments and in particular of asymmetries can be obtained from the 
correlator ratios i?'™". Indeed, the values extracted from the ratios (summarized in Table V) have smaller errors than 
those from the unconstrained analysis. The main reasons for this improvement are that we do not have to determine 
the energy E{p) and normalization constant Zn{p) for nonzero spatial momenta as both drop out in the constrained 
analysis. This reduces also the statistical noise as the nucleon correlator with smeared source and sink is not involved 
anymore in the data analysis. 

The normalization constants and Ai in Table V were determined by performing a joint fit of all relevant 
correlators. This approach is equivalent to the FC analysis method. However, as the correlators with higher momenta 
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are not involved the obtained results have smaller errors compared to the FC analaysis. Our values for a = —0.0091 ± 
0.0002st±0.0003sys and /3 = 0.0090±0.0002st±0.0003sys obtained from A,/mAr at /3 = 5.40 (see Table V) are consistent 
within the errors with the recent results a = -0.0112 ± 0.0012st ± 0.0022sys and /5 = 0.00120 ± O.OOlS^t ± 0.0023sys 
from simulations with 2+1 flavors of domain- wall fermions [46]. 





FIG. 7: The effect of different chiral extrapolations is demonstrated in the case of i?'^"" where in (a) a linear fit is performed 
and in (b) a quadratic one. In the lower plots we show the chiral extrapolation of the asymmetry 7?^"° — i?"^" (c) and the sum 
i?^'"' + F^^'^ (d). All the plots contain also one, two and three sigma error bands of the corresponding fits. 




FIG. 8: Linear chiral extrapolation of the second moment ratio J?^^" (a) and of the asymmetry J?^*^" — ii"^" (b) as obtained 
from the constrained analysis with one, two and three sigma error bands of the corresponding fits. 

In principle one can calculate similar ratios for correlators involving 
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instead of using 

However, this leads to statistical errors which are about three times larger. 

To illustrate the dependence of i?^"" on the pion mass we present in Figs. 7(a) and 7(b) linear and quadratic chiral 
extrapolations of this quantity. As i?"^'^ exhibits a similar behavior, but with opposite slope, the deviation from a 
linear dependence is amplified in the asymmetry — R°^° (Fig- 7(c)). On the other hand, this leads to linear 
behavior of + (Fig. 7(d)). Thus, due to momentum conservation one expects also linear behavior for R^'~'^, 
which is indeed observed in our data. Of course deviations from linear behavior are also possible for all other moments. 
However, they seem to be smaller than present statistical errors. Comparing the chiral extrapolations in Fig. 5(a) to 
Fig. 7(c) and in Fig. 5(b) to Fig. 7(a) reveals the increased accuracy of the constrained analysis. 

This increase of accuracy is even more important for higher moments. From Figs. 8(a) and 8(b) it is obvious that 
the improvement for the second moments allows us not only to determine the moments but also the more interesting 
asymmetries. Even more, with the help of the constraints (24) the moments f^^"", 0"^°, ^^^•^ can be calculated from 
the other second moments and the first moments. Our results are fully consistent with the direct determination. This 
approach can be particularly advantageous in the calculation of the third moments as one can then dispense with the 
evaluation of 





/3 


5.40 


5.29 


/jv/rn^ • 10^ 


3.573(69)(33)(61) 


3.392(68)(178)(74) 


-Ai/mjv 


• 10='[GcV] 41.29(74)(45)(113) 


42.32(81)(277)(133) 


A2/miv 


• 10^[GcV] 


81.27(149)(90)(221) 


83.90(167)(599)(261) 




100 


0.3638(11)(68)(3) 


0.3549(11)(61)(2) 




= ,^010 


0.3023(10)(42)(5) 


0.3100(10)(73)(1) 




0.3339(9)(26)(2) 


0.3351(9)(11)(2) 


^"><> 


-0001 


0.0300(23)(93)(1) 


0.0199(23)(46)(4) 




- 0010 


0.0313(17)(12)(7) 


0.0251(16)(84)(2) 




,011 


0.0724(18)(82)(70) 


0.0863(23)(97)(74) 


,/,ioi 


= ^101 


0.1136(17)(32)(21) 


0.1135(23)(3)(33) 


^110* 


0.0937(16)(3)(38) 


0.0953(21)(58)(31) 




200 


0.1629(28)(7)(68) 


0.1508(38)(213)(64) 


^020* 


= ipO^O 


0.1289(27)(37)(51) 


0.1207(32)(43)(56) 




,002 


0.1488(32)(77)(73) 


0.1385(36)(47)(64) 


0110 


- 0011 


0.0211(27)(78)(32) 


0.0075(33)(69)(44) 


0101 


-0110 


0.0204(21)(134)(50) 


0.0172(29)(82)(57) 


^200 


- 0020 


0.0321(33)(69)(55) 


0.0335(43)(26)(78) 


0002 


-0020 


0.0193(24)(32)(42) 


0.0170(36)(8)(56) 



TABLE V: The results for (j^'™" and the relevant asymmetries as obtained from the chirally extrapolated ratios i^'™" in the 
MS renormalization scheme at 4GeV^. The values marked by a star were used in the analysis of the corresponding asymmetries 
to determine the overall normalization. The first error is the combined statistical error of the moments and renormalization 
matrices dominated by the statistical uncertainties of the moments. The second (third) errors are the systematic uncertainties 
due to the chiral extrapolation (renormalization). 

Our data do not allow us to perform a continuum extrapolation. However, the fact that the /5 = 5.29 and /3 = 5.40 
results are compatible with each other indicates that its effect would be small. Thus we take the data from our finer 
lattice (/? = 5.40) as our final numbers. For convenience we summarize in Table VI the corresponding moments (^'™" 
at two different renormalization scales as obtained from the P — 5.40 results in Table V. The change of scales has 
been performed in the one-loop approximation with Aj^ — 226 MeV. For this purpose the moments (/j'™", being 
not multiplicatively renormalizable, had to be expressed as linear combinations of quantities that are multiplicatively 
renormalizable, at least in the one-loop approximation, i.e., the coefficients c„/ to be introduced in the next section. 
Their values (and hence also the values of the moments c^'"*" at the new scale) depend somewhat on the set of moments 
ii^g^i gj.Q used as an input. We employed here the set 1 of moments defined in the following section. 

V. MODELLING THE NUCLEON DISTRIBUTION AMPLITUDE 

Since the available nonperturbative information on the nucleon DA comes in the form of a few first moments, 
it is tempting to choose a model which is polynomial in momentum fractions at the reference scale hq. A natural 
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Asymptotic 


^2^4 GoV^ 


= 1 GcV^ 


/„ • 102[GoV"] 
-Ai • 10^[GcV^] 
A2 • 10^[GoV^] 




3.144(61)(83) 
38.72(76)(148) 
76.23(139)(291) 


3.234(63)(86) 
35.57(65)(136) 
70.02(128)(268) 




i ^ 0.333 
i « 0.333 
i « 0.333 


0.3936(34)(126) 

0.3023(10)(47) 

0.3041(29)(96) 


0.3999(37)(139) 

0.2986(11)(52) 

0.3015(32)(106) 


^011 


i « 0.143 
i « 0.143 
i « 0.143 
^ « 0.095 
^ « 0.095 
^ « 0.095 


0.1788(53)(179) 

0.1289(27)(88) 

0.1310(95)(324) 

0.0659(74)(266) 

0.1072(35)(128) 

0.1076(56)(182) 


0.1816(64)(212) 

0.1281(32)(106) 

0.1311(113)(382) 

0.0613(89)(319) 

0.1091(41)(152) 

0.1092(67)(219) 



TABLE VI: Moments vs'"" as obtained from the independent subset 0°^^ 0^°" and at /? = 5.40 in Table V at 

two different scales ij? = 4 GeV^ and /j,^ = 1 GeV^ in the MS renormalization scheme. 



choice corresponds to the (truncated) expansion in contributions of multipHcatively renormaHzable (to leading order) 
operators of increasing dimension [12, 48]: 

ip{Xi,^?) = 120X1X2X3 ^ ^ CnlPnl{Xi) ( —, r j . (80) 

n=0 ;=0 \'^s\^^0) J 

Here the first subscript, n — 0,...,N, is the total number of covariant derivatives in the corresponding operator 
and simultaneously the order of the polynomial P„;(xi). The second subscript, I — 0, ... ,n, enumerates independent 
local operators of the same dimension D = n + 3. In this way the scale dependence becomes particularly simple and 
the functional form is preserved under renormalization in one-loop accuracy. In addition, thanks to the conformal 
symmetry of QCD Lagrangian, the polynomials Pni{xi) are mutually orthogonal with respect to the SL{2,M.) scalar 
product 

[dx] X1X2X3 Pmk{Xi) Pnl{Xi) (X SmnSkl ■ (81) 

By this reason, the set of moments 0'™", I + m + n < 2, calculated in this work is sufficient to determine uniquely 
all coefficients in (80) up to = 2, i.e., to second order in the quark momentum fractions. Contributions of higher 
order polynomials correspond to higher dimension operators and can be added when the corresponding information 
becomes available. 

In the literature there seems to be no standard convention for the normalization of the polynomials P„;(xi) so we 
choose the simplest expressions (cf. [12, 48]): 

V5(xi,X2,X3,^ ) =120xiX2X3|l + cio(/io)(a;i- 2x2 + X3)L''Po 

20 p „ 14 

+ cii(/xo)(xi - X3)L3ft> + c2o(mo) [1 + 7(x2 - 2x1X3 - 2xl)\ L-'^o 

40 r in 32 'l 

+ C2i(Aio) (1 - 4x2) {xi - X3) L«'^o +C22(/io) [3 - 9x2 + 8x2 _ 12x1x3] 1 , (82) 

where 

L^^4^^, Po^n-lnp. (83) 
The scale dependence of the normalization constant is to this accuracy 

fNif-i) ^ fN{^^o)L^" ■ (84) 



18 



given 


in terms of the moments (^''"" 


as 






ClO = 






0°"i) - 2) , 








Cll = 


63 , 












C20 = 


-—(0200 


+ 0O"2 + 3^101) _^ 
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(4- 


1- cioj 


C21 = 


378 












C22 = 


126 


(20200 _ 


f 20°"2^0i"i) - 


21 

y 


(4 4 





Note that for N = 2 there are five independent coefficients c„;, which is also the number of independent 
^imn Jqj, I ji < 2 due to the constraints (24). In the above expressions we have chosen 0"'^™, 0°°"'^, < 

and cfP^"^ to be the independent subset. 

O.OaI.O 



(85) 

(86) 

(87) 
(88) 
(89) 



moments 

aoi ^200 




1.0 II 0.0 

0.0 0.2 0.4 X J 0.6 0.8 1.0 

(a) 



1.0 0.0 
0.0 0.2 0.4 X J 0.6 0.8 1.0 

(b) 




0.0 0.2 0.4 JC, 0.6 0.8 1.0 

(c) 



(d) 



FIG. 9: Barycentric contour plot of the leading-twist distribution amplitude ip{xi,X2,X3, fj,'^) in the limit oi — fi^ ^ 00 (a) 
and at fi^ — 4 GeV^ (b-d) using expansion (80) as obtained from the f3 = 5.40 moments presented in Table V. The asymmetry 
caused by the first moments only (A^ = 1) is illustrated in (b), while in (c-d) we took into account also the second moments 
(A'' = 2). In (c) we have used set 1 and in (d) set 2 as described in the text. The lines of constant xi, X2 and X3 are parallel to 
the sides of the triangle labelled by X2, X3 and xi, respectively. 



Our final results for the coefficients c„/ at the renormalization scale fi^ — 4GeV2 as obtained from the (3 — 5.40 
moments presented in Table V are collected in Table VII. As the central values for the moments 0'™" with l+m+n = 2 
do not fulfill the constraint (24) exactly, the values of C20, C21, C22 depend on the set of moments 0'™" that are used 
as an input. To illustrate this effect, we show two sets of the coefficients obtained from 0^°^, fli^oo^ ,^002 ^gg^. ^^^^ 
0^"^, 0°^^, 0^^" (set 2). The difference between the two sets is, of course, part of the uncertainty of the calculation. 
We estimate the overall uncertainty to be about 30% for cio, Cn, of order 50% for C20, C21 and a factor of two for 022- 
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Set 1 


Set 2 


ClO 
Cll 


0.326 
0.940 


0.326 
0.940 


C20 
C21 
C22 


-0.872 
-3.130 
0.405 


-0.687 
-5.210 
0.036 



TABLE VII: Central values of the coefficients Cni in the expansion (82) at the renormalization scale fi^ — A GeV^ as obtained 
from the /3 = 5.40 moments presented in Table V. 



X2 

0.5 0.4 0.3 0.2 0.1 




0.1 0.2 0.3 0.4 0.5 



FIG. 10: The model function ^p{xi) for xi at xg — 0.5 with its statistical uncertainty. 

The resulting shape of the nucleon DA is illustrated in Fig. 9. The asymptotic DA corresponding to the leading term 
in the expansion (82) is shown in Fig. 9(a). It is totally symmetric in the three quark momentum fractions. The model 
obtained by adding the terms proportional to cio and cn is presented in Fig. 9(b). Compared to the asymptotic case, 
the maximum is shifted towards larger values of xi indicating that the first quark carries a larger fraction of the proton 
momentum. Finally, for the plots in Figs. 9(c) and 9(d) we add contributions of the second order polynomials {n = 2), 
using the coefficients C20, C21, C22 from the first and the second set in Table VII, respectively. The difference is in fact 
not too large and the effect is the same in both cases: The maximum is smeared out forming two local maxima and one 
local minimum. While the model function from set 2 exhibits an approximate symmetry ip{xi, X2, X3) « '^{xi, X3, X2), 
this property is less obvious in the case of set 1. However, the general pattern is preserved. 

Whether the change in the shape of the DA caused by adding the second-order polynomials is of phenomenological 
significance can only be investigated in a dedicated study, which goes beyond the scope of this work. Note, however, 
that in full analogy to usual quantum mechanics, the quality of an approximation to the wave function has to be 
measured with respect to the scalar product of the appropriate Hilbert space, in our case Eq. (81), and not pointwise 
in, e.g., the momentum fraction representation. 

In Fig. 10 we show ip{xi) at x^ = 0.5 as a function of xi {12) together with the statistical error in order to give an 
impression of the corresponding uncertainty. The effect of choosing different subsets of (^'™" with I + m + n = 2 is 
demonstrated in Fig. 11 where we plot the difference of ip{xi) for set 1 and set 2. 

VI. SUMMARY AND CONCLUSIONS 

We have evaluated the first few moments of the leading-twist nucleon DA in lattice QCD. Along with these moments 
we have determined the nucleon couplings to local subleading (twist-four) operators. The required correlators have 
been computed on gauge field configurations generated by the QCDSF/DIK collaborations using two dynamical 
flavors of clover fermions. The necessary renormalization matrices have been calculated nonperturbatively, including 
the mixing with operators containing total derivatives. 

We have worked with two different gauge couplings corresponding to (3 = 5.29 and f3 — 5.40. For the lattice sizes and 
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FIG. 11: Systematic uncertainty due to the choice of the independent subsets of 0'™" with I + m + n = 2 (for details see text). 
The lines of constant xi, X2 and are parallel to the sides of the triangle labelled by X2, x-i and x\, respectively. 

quark masses see Table II. As our final numbers we take the data from our finer lattice (/3 = 5.40). The results for the 
moments of the leading-twist DA at two different renormalization scales are presented in Table VI. The corresponding 
coefficients in the expansion of the DA in a basis of orthogonal polynomials are given in Table VII. Truncating this 
expansion at the second order, we obtain a model of the DA which is plotted in Fig. 9. Our error estimates include 
statistic and known sources of systematic uncertainties, but still have to be considered with some caution. The largest 
uncertainty is caused by the chiral extrapolation. We expect that it will be reduced in the relatively near future when 
simulations with smaller pion masses on larger lattices become available. 

Our value for the coupling /^r, which determines the normalization of the leading-twist nucleon DA, appears to be 
approximately 40% below the corresponding QCD sum rule estimates [8, 9, 19]. If confirmed, this result would deal 
yet another blow at the hopes to calculate the nucleon magnetic form factor at realistic momentum transfers within 
perturbative QCD. At the same time, the twist-four couplings Ai and A2, which are related to the normalization of 
subleading twist-four DAs, turn out to be in agreement with other estimates. These constants are relevant, e.g., for 
the description of form factors involving a helicity fiip within perturbative QCD [2.3] and also for soft (end-point) 
corrections to the form factors in the light-cone sum rule approach [24, 49]. The same constants enter the effective 
baryon chiral Largangian and can be used to estimate the proton life time within GUT models. 

The results we have obtained for the first moments of the nucleon DA are consistent with the conventional picture 
that the valence u-quark with helicity parallel to that of the proton carries the largest fraction of its momentum, but 
the effect seems to be less pronounced compared to the corresponding QCD sum rule calculations [8, 9, 19]. Our 
numbers, however, are compatible with those extracted from the fits to the electromagnetic proton form factors within 
the light-cone sum rule approach [49]. 

Our calculation of the second moments of the DA indicates the presence of considerable second-order contributions in 
the expansion in terms of orthogonal polynomials. Qualitatively, these contributions smear out the maximum forming 
two local maxima and one local minimum (see Figs. 9(c) and 9(d)). The investigation of the phenomenological conse- 
quences of these and other features of our model DA, such as the approximate symmetry i^{xi, a;2, X3) « '^{xi, 2:3, 2:2), 
requires a dedicated study, which goes beyond the scope of the present work and will be presented elsewhere. 
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APPENDIX A: DIRAC MATRICES IN WEYL REPRESENTATION 

We have used the following representation of the Euclidean Dirac matrices: 

/ i \ /0001\ /00i0\ /OOlO 

(OOiOl (00-101 (000 -il (0001, 

T'l = ' ^'^\o -1 ' ^'^\-^ ' - K q q q 1 (Al) 
V-iOOO/ VlOOO/ \0 i J \0100 

with 

/-I \ 

75 = 71727374 = I n 0^10' ^ I ^^^^"^ ~ "^"^'''^ ' 

1/ 

The charge conjugation matrix has been chosen as 

C = 7274 . (A3) 



APPENDIX B: OPERATOR RELATIONS FOR LEADING-TWIST DISTRIBUTION AMPLITUDES 



In the following we give the relations between the operators whose matrix elements define moments of the leading- 
twist DA of spin-1 /2 baryons (DA operators) and the irreducible operators that appear in the general group-theoretical 
classification in [ 'u]. The relations are written for general quark flavors /, 5, h\ the proton case is obtained by the 
replacement f,g ^ u, h ^ d and the appropriate symmetrization to single out the contribution of isospin 1/2. 

The total symmetrization in space-time indices denoted by the curly brackets, e.g., 

y{23} ^ i_(y23_^y32) 

reflects the leading-twist projection. For example, the moment 1/"°^ is calculated from 

^e^"^ {[fmi {c^^u imtp i^DsiiMoM + ifmi (cisu imtp [^D,{^,h{o))t) . (bi) 

In the notation used below, it is not indicated explicitly on which quark the derivatives act in the operators on the 
right-hand side. However, it is always implied that the positions of the derivatives are the same on both sides of the 
equations. 



0th moment 



{BZ",~BZ",-BZ%B°^?) =4 (7374 [72ri +7ir2]) (B2) 
{BZ", -BZ", -BZ"o,BZ") -\ (7172 + 73r^]) (B3) 

{BZ", -BZ^, -BZ',BZ"i) (7172 ba' - 73r^] + 7374 [71^^ - (B4) 

The ( ySg"°) operators from the symmetry class — 1-+ (H h) are obtained from the above operators by replacing 

T on the right hand side hy V + A {V ~ A). 
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1st moments 



{B\y, B\y, -B\y) =^ (27473T{i2} + 7472T{i3} + j^^.T^'^} + 7471^^''^ + liTs^^^*}) (B5) 

-4T' ^tT' -^7™S) (27i72T{=^4> + 7472T{i3} + ^g^^T^"} + 7i74T{23} + ^.j^T^^^y^ (B6) 

(476", ^^r, 475, 4™^) =i (7274T{"} + 7273T{"} + ^i74T{23} + 7i73T{24}) (B7) 

The 65""' ( Sg'i") operators from the symmetry class £> — 1-+ {D -\ h) are obtained from the above operators by 

replacing T on the right hand side by V + ^ ( V — ^) . 

2nd moments 



{-Biy, -e^X' fi^X' ^eT) =X (^4r{123} + ^3T{124} + ^^^{134} ^ ^^^{234}^ (gg) 

The ,64™" ( operators from the symmetry class DD — 1-+ {DD -\ h) are obtained from the above operators 

by replacing T on the right hand side by V + ^ {V — A). 

APPENDIX C: RAW LATTICE RESULTS 



In this appendix we collect the results of the linear (in m^) extrapolation of our bare lattice data. The errors given 
are purely statistical. 







13 = 


5.29 




/3 = 5.40 






all 
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all 






# 


xVd.o.f 


# 


xVd.o.f 


# 


xVd.o.f 


fN/ml ■ 103 


4.088(77) 


6.563 


4.53(14) 


0.555 


4.287(74) 


0.658 


-Xi/ruN ■ 103[GeV] 


60.80(106) 


19.31 


69.28(176) 


6.209 


59.40(95) 


1.060 


-Ai • 103[GeV2] 


77.33(149) 


18.46 


82.24(209) 


3.484 


72.86(135) 


1.901 


Az/m.Y • 103[GeV] 


129.76(214) 


19.98 


141.53(360) 


4.928 


119.16(191) 


1.498 


A2 ■ 103[GeV2] 


158.00(315) 


18.31 


168.30(428) 


2.388 


146.48(270) 


2.716 




0.2987(49) 


1.125 


0.315(10) 


0.033 


0.2939(59) 


1.384 


^010 


0.2746(48) 


0.768 


0.263(11) 


0.765 


0.2719(62) 


0.335 


<^ooi 


0.2840(48) 


1.566 


0.271(11) 


2.555 


0.2740(60) 


0.972 


0011 


0.0647(37) 


0.276 


0.0633(87) 


0.711 


0.0646(44) 


1.831 


0101 


0.0606(39) 


0.821 


0.067(12) 


0.744 


0.0688(55) 


1.057 


0110 


0.0651(32) 


0.712 


0.0592(79) 


0.445 


0.0707(39) 


0.610 


0200 


0.1149(54) 


2.367 


0.146(14) 


0.597 


0.1126(68) 


5.534 


0020 


0.0922(50) 


0.717 


0.096(12) 


1.908 


0.0949(61) 


0.288 


0002 


0.1067(54) 


0.944 


0.108(13) 


2.729 


0.1060(64) 


0.114 



TABLE VIIL Linear extrapolations of FC (unconstrained) results to the physical point using all available lattice ensembles 
(all) and 24^ x 48 lattices only (24) for /3 = 5.29. The x^/d.o.f refers to the linear chiral extrapolation. 
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/3 = 5.29 = 5.40 



all 24 all 





# 


xVd.o.f 


# 


xVd.o.f 








4.396(99) 


2.417 


4.67(19) 


1.208 


4.517(96) 


0.342 


ylOO ^010 


0.308(13) 


0.416 


0.298(35) 


0.027 


0.298(19) 


0.966 


. 1 nn . ni n 

A \-\J\J 


0.0133(40) 


2.495 


0.046(13) 


0.038 


0.0196(64) 


0.960 


inn n 1 n 


0.307(12) 


0.425 


0.297(25) 


0.263 


0.300(16) 


0.483 




0.324(16) 


0.352 


0.360(49) 


0.001 


0.323(24) 


0.777 


010 010 rn 001 

¥■ = = -1 


0.286(12) 


1.636 


0.248(26) 


0.550 


0.276(17) 


0.446 


^001 _ yOOl 


0.289(15) 


1.892 


0.229(37) 


1.532 


0.280(21) 


0.399 


0100 _ 0010 


0.0194(49) 


2.230 


0.054(15) 


0.056 


0.0258(77) 


0.928 


.inn n n 1 

1 U U L 


0,0076(39) 


2.017 


0.036(14) 


1.011 


0.0129(66) 


1.291 




0.0114(41) 


0.07IJ 


O.0i(itl:!j 


i.719 


0.(J14il,(i()) 


2.118 




0.0698(56) 


0.197 


0.072(17) 


0.228 


0.0676(69) 


0.260 


-nil -ini 
^Ull _ ^^"^^ 


-0.0006(49) 


0.038 


0.000(15) 


0.004 


0.0022(60) 


1.063 




0.0689(44) 


0.395 


0.068(12) 


0.035 


0.0707(54) 


0.580 




0.0709(85) 


0.068 


0.076(27) 


0.061 


0.064(11) 


0.533 


ini ,ini iin 

<^^"i = (p'-"^ = T " 


0.0699(62) 


0.428 


0.071(18) 


0.135 


0.0673(67) 


0.504 


^110 ^ yllO 


0.0637(79) 


0.149 


0.064(24) 


0.101 


0.077(10) 


0.049 


. 1 m .mi 


0.0012(62) 


0.068 


0.006(19) 


0.023 


0.0005(73) 


1.711 




0.0025(45) 


0.048 


0.004(15) 


0.096 


-0.0042(62) 


0.246 


(/)^"^ - 0^^° 


-0.0001(47) 


0.155 


0.005(17) 


0.383 


-0.0036(62) 


0.627 


y200 _ y020 


0.1059(78) 


0.557 


0.129(22) 


0.015 


0.115(10) 


2.034 




0.0132(59) 


0.698 


0.036(18) 


0.131 


0.0195(81) 


1.812 


^^^200 _ rp020 


0.1108(79) 


0.576 


0.119(19) 


1.336 


0.1203(89) 


1.450 


200 


0.117(12) 


0.739 


0.165(37) 


0.006 


0.134(16) 


2.305 


^020 _ 0020 _ y002 


0.0913(73) 


0.261 


0.097(19) 


0.590 


0.0963(93) 


0.646 


^002 ^ ^002 
0200 _ 0020 


0.096(12) 
0.0206(68) 


0.724 
0.406 


0.066(35) 
0.039(21) 


1.320 
0.001 


0.106(15) 
0.0300(97) 


0.279 
1.864 


0200 _ 0002 
0002 _ 0020 


0.0060(61) 
0.0114(55) 


0.847 
0.291 


0.032(20) 
0.005(19) 


0.601 
0.757 


0.0092(83) 
0.0215(80) 


1.380 
0.438 



TABLE IX: Linear extrapolations of PC (unconstrained) results to the physiral point using all available lattice ensembles (all) 
and 24'^ x 48 lattices only (24) for /3 = 5.29. The x^/d.o.f refers to the linear chiral extrapolation. 







(3 = 5.29 




= 5.40 




# 


x'/d.o.f 


# 


xVd.o.f 


fN/m% ■ 10=* 


4.215(85) 


1.878 


4.395(85) 


0.267 


-Xi/mN ■ 10=*[GeV] 


51.10(117) 


10.57 


60.35(108) 


0.184 


A2/mjv • 10='[GeV] 


125.75(25) 


10.54 


120.80(216) 


0.403 




0.3286(12) 


7.559 


0.3358(11) 


6.115 


0010 ^ ^010 


0.2943(9) 


8.530 


0.2891(9) 


6.960 




0.3164(9) 


1.112 


0.3155(9) 


1.312 


0100 _ 0010 


0.0350(20) 


9.960 


0.0468(19) 


7.732 


0100 _ 0001 


0.0126(19) 


3.996 


0.0206(18) 


3.300 


0001 _ 0010 


0.0225(14) 


3.315 


0.0263(14) 


2.526 


</,°" 


0.1113(26) 


3.593 


0.0932(19) 


1.544 




0.1148(26) 


0.370 


0.1124(18) 


0.287 




0.1085(22) 


1.716 


0.1034(16) 


0.135 


0200 


0.1820(44) 


4.176 


0.1924(30) 


0.338 


0020 ^^020 (^j 


0.1489(35) 


0.363 


0.1539(28) 


0.265 


0002 


0.1728(42) 


1.677 


0.1801(36) 


0.856 




0.0042(39) 


2.489 


0.0200(27) 


0.900 


0110 _ 0011 


0.0042(34) 


0.636 


0.0100(25) 


0.775 


0101 - 0110 


0.0053(29) 


1.159 


0.0094(20) 


0.257 


0200 _ 0020 


0.0367(48) 


1.515 


0.0364(35) 


0.514 


0200 _ 0002 


0.0076(59) 


1.763 


0.0115(39) 


0.810 


0002 _ 0020 


0.0230(39) 


1.010 


0.0255(24) 


0.597 



TABLE X: Linear extrapolations of (^'"'" and asymmetries to the physical point as obtained from the constrained analysis 
using all available ensembles. The x^/d.o.f refers to the linear chiral extrapolation. The values denoted by the ★ were used to 
determine the absolute normalization of the associated asymmetries. 
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